########################
# Manipulating Remittances
# Christopher A. Culver
# Last updated: 25 May 2021
########################

setwd("C:/Users/Christopher.Culver/OneDrive - afacademy.af.edu/Manipulating Remittances Submission")

install.packages("foreign")
install.packages("countrycode")
install.packages("reshape")
install.packages("DataCombine")
install.packages("interplot")
install.packages("survival")
install.packages("TTR")
install.packages("VIM")
install.packages("mice")
install.packages("plm")
install.packages("stargazer")
install.packages("car")
install.packages("tidyr")
install.packages("dynpanel")
install.packages("ecm")
install.packages("multcomp")
#install.packages("VIMGUI")
install.packages("psych")
install.packages("sjPlot")
install.packages("margins")
install.packages("lme4")

library(foreign)
library(countrycode)
library(reshape)
library(DataCombine)
library(interplot)
library(survival)
library(TTR)
library(VIM)
library(mice)
library(plm)
library(stargazer)
library(car)
library(tidyr)
library(dynpanel)
library(ecm)
library(multcomp)
#library(VIMGUI)
library(psych)
library(sjPlot)
library(margins)
library(lme4)

rm(list = ls())

##########################################

POLITY <- read.csv("p4v2015.csv") #up to 2015
GWF <- read.csv("GWFtscs2015.csv") #up to 2015
WDI <- read.csv("WDI.csv") #up to 2015
PRIO <- read.csv("PRIOv4.csv") #up to 2014
OIL <- read.csv("Ross Oil Data.csv") #up to 2006
KA <- read.csv("kaopen_2013.csv") #up to 2013
REER <- read.csv("REER_Annual_67.csv") #up to 2015
ICTD <- read.csv("ICTDGRD_June2016_CentralGeneralMerged.csv") #1980-2013
FLIGHT <- read.csv("Captial Flight Data.csv")  
CBI <- read.dta("CBI dataset Garriga_small.dta")
 

##########################################

names(GWF)[names(GWF)=="cowcode"] <- "ccode"
GWF <- GWF[which(GWF$year>1974),]
Global <- merge(GWF, POLITY , by=c("ccode", "year"), all.x=T)

names(WDI)[names(WDI)=="�..Country.Code"] <- "wbcode"
names(WDI)[names(WDI)=="Time.Code"] <- "year"
names(WDI)[names(WDI)=="BX.TRF.PWKR.DT.GD.ZS"] <- "remitbygdp"
names(WDI)[names(WDI)=="BX.TRF.PWKR.CD.DT"] <- "remit"
names(WDI)[names(WDI)=="NY.GDP.PCAP.CD"] <- "gdppc"
names(WDI)[names(WDI)=="NY.GDP.MKTP.KD.ZG"] <- "gdpgrowth"
names(WDI)[names(WDI)=="SP.POP.TOTL"] <- "pop"
names(WDI)[names(WDI)=="SM.POP.NETM"] <- "migration"
names(WDI)[names(WDI)=="PX.REX.REER"] <- "reer2"
names(WDI)[names(WDI)=="BX.KLT.DINV.WD.GD.ZS"] <- "fdibygdp"
names(WDI)[names(WDI)=="BN.KLT.DINV.CD"] <- "fdi"
names(WDI)[names(WDI)=="DT.ODA.ALLD.KD"] <- "aid"
names(WDI)[names(WDI)=="FP.CPI.TOTL.ZG"] <- "cpi"
names(WDI)[names(WDI)=="NY.GDP.DEFL.KD.ZG"] <- "inflation"
names(WDI)[names(WDI)=="NE.EXP.GNFS.ZS"] <- "exports"
names(WDI)[names(WDI)=="NE.IMP.GNFS.ZS"] <- "imports"
names(WDI)[names(WDI)=="NE.TRD.GNFS.ZS"] <- "trade"
names(WDI)[names(WDI)=="FM.LBL.BMNY.CN"] <- "broadmoney"
names(WDI)[names(WDI)=="FM.AST.CGOV.ZG.M3"] <- "claimsongov"
names(WDI)[names(WDI)=="FM.LBL.BMNY.IR.ZS"] <- "broadtoreserves"
names(WDI)[names(WDI)=="FM.LBL.BMNY.ZG"] <- "broadmoneygrowth"
names(WDI)[names(WDI)=="FM.LBL.BMNY.GD.ZS"] <- "broadmoneybygdp"
names(WDI)[names(WDI)=="FI.RES.TOTL.DT.ZS"] <- "reservesbydebt"
names(WDI)[names(WDI)=="FI.RES.TOTL.MO"] <- "resevemonths"
names(WDI)[names(WDI)=="BN.CAB.XOKA.GD.ZS"] <- "currentacct"
names(WDI)[names(WDI)=="NV.IND.MANF.ZS"] <- "manufacturing"
names(WDI)[names(WDI)=="NY.GDP.MKTP.CD"] <- "gdp"
names(WDI)[names(WDI)=="PA.NUS.PPPC.RF"] <- "nationalprice"

WDI$year <- gsub("YR", "", WDI$year)
WDI$year <- as.numeric(WDI$year)
WDI <- WDI[which(WDI$year !=""),]
WDI$ccode <- countrycode(WDI$wbcode,"wb","cown",warn=T)
for(i in 1:length(WDI$ccode)){
 if(WDI$wbcode[i]=="YEM"){
   WDI$ccode[i] <- 678
 }
}
Global <- merge(WDI, Global, by=c("ccode", "year"))
#Congo GDP data in WDI listed as 0 instead of NA#
Global[which(Global$ccode==490 & Global$year < 1991),"gdppc"] <- NA 


#########
names(REER)[names(REER)=="Updated..25.November.2015"] <- "year"
names(REER)<- gsub("REER_67_", "", names(REER))
REER <- melt(REER, id="year")
REER <- REER[which(REER$year != "NA"),]
names(REER)[names(REER)=="variable"] <- "iso2"
names(REER)[names(REER)=="value"] <- "reer"
REER$ccode <- countrycode(REER$iso2,"iso2c","cown",warn=T)
Global <- merge(Global, REER, by=c("ccode", "year"), all.x=T)
Global$reer <- as.numeric(Global$reer)

OIL <- OIL[,c("cty","id","year","oil_gas_rentPOP")]
OIL$ccode <- countrycode(OIL$id,"wb","cown",warn=T)
Global <- merge(Global, OIL, by=c("ccode", "year"), all.x=T)

names(KA)[names(KA)=="ccode"] <- "wbcode"
KA$cn <- NULL
Global <- merge(Global, KA, by=c("wbcode", "year"), all.x=T)

names(ICTD)[names(ICTD)=="Calendar.year..nearest."] <- "year"
names(ICTD)[names(ICTD)=="ISO"] <- "wbcode"
names(ICTD)[names(ICTD)=="Taxes.on.goods.and.services..of.which.Taxes.on.Sales"] <- "salestax"
names(ICTD)[names(ICTD)=="Taxes.on.international.trade.and.transactions.Of.which.Export"] <- "exporttax"
names(ICTD)[names(ICTD)=="Other.taxes"] <- "othertax"
names(ICTD)[names(ICTD)=="Non.Resource.Component.of.Non.Tax"] <- "otherrevenue"
names(ICTD)[names(ICTD)=="Revenue.excluding.Grants..including.social.contributions."] <- "totalrevenue"
Global <- merge(Global, ICTD, by=c("wbcode", "year"), all.x=T)

Global$migration <- as.numeric(Global$migration)
cow.codes <- unique(Global$ccode)
 for(i in 1:length(cow.codes)) {
   subset <- Global[which(Global$ccode==cow.codes[i]),c("migration","ccode","year")]
   subset <- FillDown(subset,'migration')
   Global[which(Global$ccode==cow.codes[i]),"migration"] <- subset$migration
 }

FLIGHT <- melt(FLIGHT, id=c("ccode","Country"),variable_name = "year") #doesn't work like it used to, code doesn't make sense but works
names(FLIGHT)[names(FLIGHT)=="variable"] <- "year"
FLIGHT$year <- gsub("X", "", FLIGHT$year)
FLIGHT$year <- as.numeric(FLIGHT$year)
FLIGHT <- FLIGHT[order(FLIGHT$ccode, FLIGHT$year),]
Global <- merge(Global, FLIGHT, by=c("ccode", "year"), all.x=T)

CFA <- c(434,439,437,436,433,461,471,482,483,484,481,432)
Global$cfa <- 0
for(i in 1:length(CFA)) {
 for(j in 1:length(Global$cfa)){
  if(Global$ccode[j]==CFA[i]){
   Global$cfa[j] <- 1
  }
 }
}

###Mali out of CFA 1962-1983###
for(i in 1:length(Global$cfa)){
 if(Global$ccode[i]==432 & Global$year[i]>1961 & Global$year[i]<1984){
   Global$cfa[i] <- 0
 }
###Eq Guinea in CFA after 1984 ###
 if(Global$ccode[i]==411 & Global$year[i]>1984){
   Global$cfa[i] <- 1
 }
###Guinea-Bissau in CFA after 1996 ###
 if(Global$ccode[i]==404 & Global$year[i]>1996){
   Global$cfa[i] <- 1
 }
###Madagascar in CFA before 1973 ###
 if(Global$ccode[i]==580 & Global$year[i]<1973){
   Global$cfa[i] <- 1
 }
###Mauritania in CFA before 1973 ###
 if(Global$ccode[i]==435 & Global$year[i]<1973){
   Global$cfa[i] <- 1
 }
}

Global$dollarized <- 0
for(i in 1:length(Global$dollarized)){
###Zimbabwe abandonded currency 2009###
 if(Global$ccode[i]==552 & Global$year[i]>2008){
   Global$dollarized[i] <- 1
 }
###El Salvador and Ecuador dollarized 2000###
 if(Global$ccode[i]==92 & Global$year[i]>1999){
   Global$dollarized[i] <- 1
 }
 if(Global$ccode[i]==130 & Global$year[i]>1999){
   Global$dollarized[i] <- 1
 }
}

#Global[,c("ccode","year","dollarized","cfa","country")]

Global$convertible <- 0
for(i in 1:length(Global$convertible )){
##China 2003##
 if(Global$ccode[i]==710 & Global$year[i]>2002){
   Global$convertible[i] <- 1
 }
##Malaysia 1997##
 if(Global$ccode[i]==820 & Global$year[i]>1996){
   Global$convertible[i] <- 1
 }
##Mexico 1997##
 if(Global$ccode[i]==70 & Global$year[i]>1996){
   Global$convertible[i] <- 1
 }
##Russia 1997##
 if(Global$ccode[i]==365 & Global$year[i]>1996){
   Global$convertible[i] <- 1
 }
##Saudi Arabia 1997##
 if(Global$ccode[i]==670 & Global$year[i]>1996){
   Global$convertible[i] <- 1
 }
##Singapore##
 if(Global$ccode[i]==830){
   Global$convertible[i] <- 1
 }
##South Africa##
 if(Global$ccode[i]==560){
   Global$convertible[i] <- 1
 }
##Spain##
 if(Global$ccode[i]==230){
   Global$convertible[i] <- 1
 }
##Thailand 1997##
 if(Global$ccode[i]==800 & Global$year[i]>1996){
   Global$convertible[i] <- 1
 }
}

#Global[,c("ccode","year","gwf_country","convertible","cfa","dollarized")]

Global$overvalue <- 0
for(i in 1:length(Global$overvalue)){
 if(is.na(Global$reer[i])){ 
   Global$overvalue[i] <- NA
 } else {
    if(Global$reer[i]>100){
      Global$overvalue[i] <- 1
    }
 }
}

#Global[,c("reer","overvalue")]

Global$controller <- 0
for(i in 1:length(Global$controller)){
 if(Global$cfa[i]==0 & Global$dollarized[i]==0 & Global$convertible[i]==0){
   Global$controller[i] <- 1
 }
}

Global$overvalued_controller <- 0
for(i in 1:length(Global$controller)){
 if(is.na(Global$overvalue[i])){ 
   Global$overvalued_controller[i] <- NA
 } else {
     if(Global$cfa[i]==0 & Global$dollarized[i]==0 & Global$convertible[i]==0 & Global$overvalue[i]==1){
     Global$overvalued_controller[i] <- 1
    }
 }
}


#Global[,c("ccode","year","cfa","convertible","dollarized","overvalue","controller","overvalued_controller","country")]

names(CBI)[names(CBI)=="cowcode"] <- "ccode"
Global <- merge(Global, CBI, by=c("ccode", "year"), all.x=T)


#####Short-cut########
Data <- read.csv("Global.csv")
Global$conflict <- Data$conflict

####Long-cut############
PRIO <- subset(PRIO, TypeOfConflict>2, select=c("Year","Location","IntensityLevel","TypeOfConflict", "GWNoA"))
Global$conflict <- 0
for(i in 1:length(Global$conflict)) {
   for(j in 1:nrow(PRIO)) {
     if(PRIO$Year[j]==Global$year[i] & PRIO$GWNoA[j]==Global$ccode[i] & Global$conflict[i]<2){
       Global$conflict[i] <- PRIO$IntensityLevel[j]
     }
   }
}
#####
write.csv(Global, "C:/Users/Christopher.Culver/OneDrive - afacademy.af.edu/Manipulating Remittances Submission/Global.csv")


##############################
rm(list = ls())
Global <- read.csv("Global.csv")

Data <- Global[,c("ccode","year","remit","gdppc","gdpgrowth","pop","remitbygdp",
  "reservesbydebt","resevemonths","currentacct","manufacturing","gdp",
  "migration","fdibygdp","fdi","aid","cpi","inflation","exports","imports","reer",
  "trade","gwf_country","polity2","gwf_duration","gwf_fail","gwf_party","gwf_fail_type",
  "reer2","oil_gas_rentPOP","kaopen","conflict","broadmoney","broadmoneygrowth",
  "salestax","totalrevenue","otherrevenue","exporttax","othertax","convertible",
  "overvalue","controller","overvalued_controller","nationalprice","lvaw_garriga")]

###removed due to error "flight"

write.csv(Data, "C:/Users/Christopher.Culver/OneDrive - afacademy.af.edu/Manipulating Remittances Submission/Data.csv")


##############################





##############################
##############################
##############################
##############################

rm(list = ls())
Data <- read.csv("Data.csv")
#nrow(Data)

names(Data)[names(Data)=="gwf_duration"] <- "time"
names(Data)[names(Data)=="gwf_fail"] <- "fail"
names(Data)[names(Data)=="gwf_party"] <- "party"
names(Data)[names(Data)=="oil_gas_rentPOP"] <- "oilraw"
names(Data)[names(Data)=="lvaw_garriga"] <- "cbi"
names(Data)[names(Data)=="remit"] <- "remitraw"
names(Data)[names(Data)=="pop"] <- "popraw"
names(Data)[names(Data)=="aid"] <- "aidraw"

Data$manufacturingtotal <- Data$manufacturing/100*Data$gdp
Data$exportstotal <- Data$exports/100*Data$gdp

Data$coup <- 0
for(i in 1:length(Data$coup)){
 if(Data$fail[i]==1 & Data$gwf_fail_type[i]!=4 & Data$gwf_fail_type[i]!=7){
   Data$coup[i] <- 1
 }
}
#Data[,c("fail","gwf_fail_type","coup")]

Data$time2 <- Data$time^2
Data$time3 <- Data$time^3
Data$trend <- (Data$year - 1974)
Data$trend2 <- Data$trend^2

Data$reer3 <- NA
Data$reer4 <- NA

PRICE <- read.csv("national price level.csv")
PRICE$ccode <- countrycode(PRICE$Country.Code,"wb","cown",warn=T)

cow.codes <- unique(Data$ccode)
for(i in 1:length(cow.codes)) {
 if(cow.codes[i]!=360 & cow.codes[i]!=490 & cow.codes[i]!=678){
  x07 <- PRICE[which(PRICE$ccode==cow.codes[i]),"X2007"]
  x11 <- PRICE[which(PRICE$ccode==cow.codes[i]),"X2011"]
  Data[which(Data$ccode==cow.codes[i]),"reer3"] <- Data[which(Data$ccode==cow.codes[i]),"reer"]*x07
  Data[which(Data$ccode==cow.codes[i]),"reer4"] <- Data[which(Data$ccode==cow.codes[i]),"reer"]*x11
 }
}

Data$reer3 <- log(Data$reer3)
Data$reer4 <- log(Data$reer4)


##############Transformations##########
Data$remit <- NA
for(i in 1:nrow(Data)) {
 if(!is.na(Data$remitraw[i]) & Data$remitraw[i]==0){
   Data$remit[i] <- 0
 } else if(!is.na(Data$remitraw[i]) & Data$remitraw[i]>0) {
   Data$remit[i] <- log(Data$remitraw[i]/Data$popraw[i])
   }
}
Data$reer <- log(Data$reer)
Data$reer2 <- log(Data$reer2)
Data$gdppc <- log(Data$gdppc)
Data$gdp <- log(Data$gdp)
Data$aid <- Data$aidraw/Data$popraw
Data$migration <- Data$migration/Data$popraw
Data$pop <- log(Data$popraw)
Data$exports <- log(Data$exports)

Data$seigniorage <- Data$broadmoneygrowth/Data$inflation

Data$oil <- NA
for(i in 1:length(Data$oil)) {
 if(!is.na(Data$oilraw[i]) & Data$oilraw[i] > 0){
   Data$oil[i] <- log(Data$oilraw[i])
 }
}

Data$salestax <- as.numeric(gsub("%", "", Data$salestax))
for(i in 1:length(Data$salestax)) {
 if(!is.na(Data$salestax[i]) & Data$salestax[i] > 0){
   Data$salestax[i] <- log(Data$salestax[i])
 }
}

Data$otherrevenue <- as.numeric(gsub("%", "", Data$otherrevenue))
for(i in 1:length(Data$otherrevenue)) {
 if(!is.na(Data$otherrevenue[i]) & Data$otherrevenue[i] > 0){
   Data$otherrevenue[i] <- log(Data$otherrevenue[i])
 }
}

for(i in 1:length(Data$seigniorage)) {
 if(!is.na(Data$seigniorage[i]) & Data$seigniorage[i] > 0){
   Data$seigniorage[i] <- log(Data$seigniorage[i])
 }
 if(!is.na(Data$seigniorage[i]) & Data$seigniorage[i] < 0){
   Data$seigniorage[i] <- -1*log(-1*Data$seigniorage[i])
 }
}

Data$aidlog <- NA
for(i in 1:length(Data$aid)) {
 if(!is.na(Data$aid[i]) & Data$aid[i] > 0){
   Data$aidlog[i] <- log(Data$aid[i])
 }
 if(!is.na(Data$aid[i]) & Data$aid[i] < 0){
   Data$aidlog[i] <- -1*log(-1*Data$aid[i])
 }
}


Data$totalrevenue <- as.numeric(gsub("%", "", Data$totalrevenue))
for(i in 1:length(Data$totalrevenue)) {
 if(!is.na(Data$totalrevenue[i]) & Data$totalrevenue[i] > 0){
   Data$totalrevenue[i] <- log(Data$totalrevenue[i])
 }
}

Data$othertax <- as.numeric(gsub("%", "", Data$othertax))
Data$otherrevenue <- as.numeric(gsub("%", "", Data$otherrevenue))
Data$other <- NA
A <- NULL
B <- NULL
for(i in 1:length(Data$othertax)){
 if(is.na(Data$othertax[i])){
  A[i] <- 0
 } else {
  A[i] <- Data$othertax[i]
 } 
 if(is.na(Data$otherrevenue[i])){
  B[i] <- 0
 } else {
  B[i] <- Data$otherrevenue[i]
 } 
} 
Data$other <- A+B
for(i in 1:length(Data$other)){
 if(Data$other[i]==0){
  Data$other[i] <- NA
 }
}


##############/Transformations##########


#############Missingness/Imputation##########
SurvivalData <- Data[,c("ccode","year","reer","remit","party","gdppc","gdpgrowth","pop","conflict","remitbygdp",
   "aid","migration","kaopen","oil")]
#aggr_plot <- aggr(SurvivalData, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(SurvivalData), 
#   cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))
#sapply(SurvivalData, function(x) sum(is.na(x)))
#names(SurvivalData)

set.seed(156)
init = mice(SurvivalData, maxit=0) 
meth = init$method
predM = init$predictorMatrix
predM[, c("ccode")]=0
meth[c("ccode","year","reer","party","pop","conflict","aid","migration","kaopen","oil")]=""
#meth[c("remit","gdppc","gdpgrowth","remitbygdp")]="norm" 
imputed = mice(SurvivalData, method=meth, predictorMatrix=predM, m=5) 
Imputed <- mice::complete(imputed)

Data$remitimp <- Imputed$remit 
Data$remitbygdpimp <- Imputed$remitbygdp
Data$gdppcimp <- Imputed$gdppc 
Data$gdpgrowthimp <- Imputed$gdpgrowth

#Data[10:20,c("ccode","year","remit","remitimp","gdppc","gdppcimp","gdpgrowth","gdpgrowthimp")]
#sapply(Data, function(x) sum(is.na(x)))

#################/Imputation#########

################LagVariables##########

Data <- Data[order(Data$ccode, Data$year),]

Data <- slide(Data, Var = "remit", TimeVar = "year", GroupVar = "ccode", NewVar = "remitlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "remitimp", TimeVar = "year", GroupVar = "ccode", NewVar = "remitimplag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "reer", TimeVar = "year", GroupVar = "ccode", NewVar = "reerlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "remitbygdp", TimeVar = "year", GroupVar = "ccode", NewVar = "remitbygdplag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "remitbygdpimp", TimeVar = "year", GroupVar = "ccode", NewVar = "remitbygdpimplag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "reer2", TimeVar = "year", GroupVar = "ccode", NewVar = "reer2lag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "party", TimeVar = "year", GroupVar = "ccode", NewVar = "partylag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "gdpgrowth", TimeVar = "year", GroupVar = "ccode", NewVar = "gdpgrowthlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "gdppc", TimeVar = "year", GroupVar = "ccode", NewVar = "gdppclag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "pop", TimeVar = "year", GroupVar = "ccode", NewVar = "poplag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "conflict", TimeVar = "year", GroupVar = "ccode", NewVar = "conflictlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "aid", TimeVar = "year", GroupVar = "ccode", NewVar = "aidlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "migration", TimeVar = "year", GroupVar = "ccode", NewVar = "migrationlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "kaopen", TimeVar = "year", GroupVar = "ccode", NewVar = "kaopenlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "oil", TimeVar = "year", GroupVar = "ccode", NewVar = "oillag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "exports", TimeVar = "year", GroupVar = "ccode", NewVar = "exportslag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "inflation", TimeVar = "year", GroupVar = "ccode", NewVar = "inflationlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "gdp", TimeVar = "year", GroupVar = "ccode", NewVar = "gdplag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "resevemonths", TimeVar = "year", GroupVar = "ccode", NewVar = "resevemonthslag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "currentacct", TimeVar = "year", GroupVar = "ccode", NewVar = "currentacctlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "polity2", TimeVar = "year", GroupVar = "ccode", NewVar = "polity2lag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "manufacturing", TimeVar = "year", GroupVar = "ccode", NewVar = "manufacturinglag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "seigniorage", TimeVar = "year", GroupVar = "ccode", NewVar = "seignioragelag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "cbi", TimeVar = "year", GroupVar = "ccode", NewVar = "cbilag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "salestax", TimeVar = "year", GroupVar = "ccode", NewVar = "salestaxlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "aidlog", TimeVar = "year", GroupVar = "ccode", NewVar = "aidloglag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "totalrevenue", TimeVar = "year", GroupVar = "ccode", NewVar = "totalrevenuelag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "othertax", TimeVar = "year", GroupVar = "ccode", NewVar = "othertaxlag", slideBy = -1, keepInvalid = T)
Data <- slide(Data, Var = "trade", TimeVar = "year", GroupVar = "ccode", NewVar = "tradelag", slideBy = -1, keepInvalid = T)


################/LagVariables##########

######################

unique(Data$year)
nrow(Data)
length(unique(Data$ccode))
summary(Data$remit)

Data$continent <- countrycode(Data$ccode, "cown", "continent", warn=T)
Data$region <- countrycode(Data$ccode, "cown", "region", warn=T)


write.csv(Data, "C:/Users/Christopher.Culver/OneDrive - afacademy.af.edu/Manipulating Remittances Submission/CompleteData.csv")
rm(list = ls())
Data <- read.csv("CompleteData.csv")


RemitData <- Data[,c("ccode","year","remit","region","continent","gdppc")]
names(RemitData)[names(RemitData)=="gdppc"] <- "income (logged gdp per capita)"
head(RemitData)


CONTROLLERS <- Data[which(Data$controller==1),]
NONCONTROLLERS <- Data[which(Data$controller==0),]
unique(NONCONTROLLERS$gwf_country)
length(unique(Data$gwf_country))

#NONCONTROLLERS[which(NONCONTROLLERS$gwf_country=="Mali"),"year"]
#Data[which(Data$gwf_country=="Niger"),"year"]
#unique(CONTROLLERS$gwf_country)

nrow(CONTROLLERS)
nrow(NONCONTROLLERS)
nrow(Data)
2914-528

plot.new()
 plot.window(xlim=c(0,200),ylim=c(0,300))
 points(Data$broadmoneygrowth,Data$inflation)
 axis(1,at=seq(0,200,20), labels=seq(0,200,20))
 axis(2,at=seq(0,300,100), labels=seq(0,300,100))
 box()
 title(ylab = "percent annual inflation",xlab = "percent broad money growth",
     main = "Figure 3\nInflation vs. Broad Money Growth")



############$SSA Case Study #########

AFRICA <- Data[which(Data$continent=="Africa"),]
SSA <- AFRICA[which(AFRICA$region!="Northern Africa"),]
SSACONTROLLERS <- SSA[which(SSA$controller==1),]
SSANONCONTROLLERS <- SSA[which(SSA$controller==0),]

SummaryDataAFRICA <- AFRICA[,c("fail","remit","remitimp","reer","seigniorage",
                       "party","gdppc","gdpgrowth","aidlog","conflict","controller")]
describe(SummaryDataAFRICA)

Y <- Surv(AFRICA$time,AFRICA$fail)
CoxModel <- coxph(Y ~ remitlag + controller + party + gdppclag + gdpgrowth + 
                    conflict + aidlag + trend + trend2,  data = AFRICA)
summary(CoxModel)
Test <- cox.zph(CoxModel)
Test
plot(Test[3])

nrow(SSA)
nrow(SSACONTROLLERS)
nrow(SSANONCONTROLLERS)
409+1045

length(unique(SSA$gwf_country))
C <- unique(SSACONTROLLERS$gwf_country)
NC <- unique(SSANONCONTROLLERS$gwf_country)
intersect(C,NC)

F <- c(mean(SSA$fail),mean(SSACONTROLLERS$fail),mean(SSANONCONTROLLERS$fail))
summary(SSA$fail)
summary(SSACONTROLLERS$fail)
summary(SSANONCONTROLLERS$fail)

d1 <- density(na.omit(SSACONTROLLERS$gdppc))
d2 <- density(na.omit(SSANONCONTROLLERS$gdppc))

d3 <- density(na.omit(SSACONTROLLERS$aidlog))
d4 <- density(na.omit(SSANONCONTROLLERS$aidlog))

d5 <- density(na.omit(SSACONTROLLERS$conflict))
d6 <- density(na.omit(SSANONCONTROLLERS$conflict))

d7 <- density(na.omit(SSACONTROLLERS$party))
d8 <- density(na.omit(SSANONCONTROLLERS$party))

par(mfrow=c(2,2),oma = c(0, 0, 3, 0))
plot.new()
plot.window(xlim=range(d1$x, d2$x),ylim=range(d1$y, d2$y))
lines(d1, col = "black", lty = 1)
lines(d2, col = "red", lty = 2)
title(main = "GDP pc (logged)", ylab = "Density")
box()
axis(1)
axis(2)

plot.new()
plot.window(xlim=range(d3$x, d4$x),ylim=range(d3$y, d4$y))
lines(d3, col = "black", lty = 1)
lines(d4, col = "red", lty = 2)
title(main = "Aid (logged)", ylab = "Density")
box()
axis(1)
axis(2)

plot.new()
plot.window(xlim=range(d5$x, d6$x),ylim=range(d5$y, d6$y))
lines(d5, col = "black", lty = 1)
lines(d6, col = "red", lty = 2)
title(main = "Conflict", ylab = "Density")
box()
axis(1)
axis(2)
legend(x = "topright",inset = 0, legend = c("Controllers", "Non-controllers"), 
       col=c("black","red"), lty=c(1,2), cex=.8)

plot.new()
plot.window(xlim=range(d7$x, d8$x),ylim=range(d7$y, d8$y))
lines(d7, col = "black", lty = 1)
lines(d8, col = "red", lty = 2)
title(main = "Single-party", ylab = "Density")
box()
axis(1)
axis(2)

mtext("Figure 4  Key variable distributions across SSA groups", outer = TRUE, cex = 1.4)

par(mfrow=c(1,1))
dev.new(width=5, height=4)
barplot(F, main="Figure 5\nSSA Regime Failure Rate", horiz=TRUE,
  names.arg=c("SSA average", "Controllers", "Non-controllers"))

names(SSACONTROLLERS)
Model1 <- clogit(fail ~ overvalue + party + gdppc + gdpgrowth + 
    conflict + aidlog +  time*party +
    time + time2 + time3 + trend + trend2 + 
    strata(ccode), data = SSACONTROLLERS)
Model1

Model2 <- clogit(fail ~ overvalue + party + gdppc + gdpgrowth + 
    conflict + aidlog +  time*party +
    time + time2 + time3 + trend + trend2 + 
    strata(ccode), data = SSANONCONTROLLERS)
Model2


length(C)
length(NC)

stargazer(Model1, Model2, 
          style="apsr", omit = c("time","trend"), dep.var.labels.include = FALSE,
          notes = "cubic time polynomials and trend term not displayed",
          model.numbers = TRUE, title = "Table 4.1: SSA Regime Failure")

SAMPLE <- SSA
SAMPLE$totaloil <- SAMPLE$oilraw*SAMPLE$popraw
SAMPLE <- SAMPLE[,c("year","remitraw","fdi","aidraw","popraw","totaloil","exportstotal","manufacturingtotal")]
#SAMPLE
SAMPLE$popraw <- as.numeric(SAMPLE$popraw)
SAMPLE$nre <- SAMPLE$exportstotal-SAMPLE$totaloil

Year <- sort(unique(SAMPLE$year))
Remit <- NA
Aid <- NA
Fdi <- NA
Pop <- NA
Oil <- NA
Exports <- NA
Manu <- NA
Nre <- NA

for(i in 1:length(Year)) {
  subset <- SAMPLE[which(SAMPLE$year==Year[i]),]
  Oil[i] <- sum(na.omit(subset$totaloil))
  Remit[i] <- sum(na.omit(subset$remitraw))
  Aid[i] <- sum(na.omit(subset$aidraw))
  Fdi[i] <- sum(na.omit(subset$fdi))
  Pop[i] <- sum(na.omit(subset$popraw))   
  Exports[i] <- sum(na.omit(subset$exportstotal))   
  Manu[i] <- sum(na.omit(subset$manufacturingtotal))  
  Nre[i] <- sum(na.omit(subset$nre))  
}

Aidpc <- Aid/Pop
Remitpc <- Remit/Pop
Fdipc <- Fdi/Pop
Oilpc <- Oil/Pop
Exportspc <- Exports/Pop
Manupc <- Manu/Pop
Nrepc <- Nre/Pop

#Aid
Aid2 <- SMA(Aid,2) #moving averages
Aid3 <- SMA(Aid,3)
Aid[41] <- NA
Aid2[41] <- NA
Aid3[41] <- NA
#cbind(Aid,Aid2,Aid3)

#Aidpc 
Aidpc2 <- SMA(Aidpc,2)
Aidpc3 <- SMA(Aidpc,3)
Aidpc[41] <- NA
Aidpc2[41] <- NA
Aidpc3[41] <- NA
#cbind(Aidpc,Aidpc2,Aidpc3)

#Remit
Remit2 <- SMA(Remit,2)
Remit3 <- SMA(Remit,3)
#cbind(Remit,Remit2,Remit3)

#Remitpc
Remitpc2 <- SMA(Remitpc,2)
Remitpc3 <- SMA(Remitpc,3)
#cbind(Remitpc,Remitpc2,Remitpc3)

#Fdi
Fdi2 <- SMA(Fdi,2)
Fdi3 <- SMA(Fdi,3)
#cbind(Fdi,Fdi2,Fdi3)

#Fdipc
Fdipc2 <- SMA(Fdipc,2)
Fdipc3 <- SMA(Fdipc,3)
#cbind(Fdipc,Fdipc2,Fdipc3)

#Oil
Oil2 <- SMA(Oil,2)
Oil3 <- SMA(Oil,3)
Oil[33:41] <- NA
Oil2[33:41] <- NA
Oil3[33:41] <- NA
#cbind(Oil,Oil2,Oil3)

#Oilpc
Oilpc2 <- SMA(Oilpc,2)
Oilpc3 <- SMA(Oilpc,3)
Oilpc[33:41] <- NA
Oilpc2[33:41] <- NA
Oilpc3[33:41] <- NA
#cbind(Oilpc,Oilpc2,Oilpc3)

#Exports
Exports2 <- SMA(Exports,2)
Exports3 <- SMA(Exports,3)
#cbind(Exports,Exports2,Exports3)

#Exportspc
Exportspc2 <- SMA(Exportspc,2)
Exportspc3 <- SMA(Exportspc,3)
#cbind(Exportspc,Exportspc2,Exportspc3)

Manu2 <- SMA(Manu,2)
Manu3 <- SMA(Manu,3)
Manupc2 <- SMA(Manupc,2)
Manupc3 <- SMA(Manupc,3)
cbind(Manu, Manu2, Manu3)

Nre2 <- SMA(Nre,2)
Nre3 <- SMA(Nre,3)
Nrepc2 <- SMA(Nrepc,2)
Nrepc3 <- SMA(Nrepc,3)
Nre[33:41] <- NA
Nre2[33:41] <- NA
Nre3[33:41] <- NA
Nrepc[33:41] <- NA
Nrepc2[33:41] <- NA
Nrepc3[33:41] <- NA
#cbind(Nrepc, Nrepc2, Nrepc3)


plot.new()
plot.window(xlim=c(1975,2015),ylim=c(-50,80))
par(mar=c(4,5,4,2.8))
lines(Year, Remitpc2, col = "red")
lines(Year, Aidpc2, lty = 3)
lines(Year, Fdipc2, lty = 2, col = "blue")
#lines(Year, Oilpc2, lty = 4, col = "brown")
#lines(Year, Nrepc, lty = 5, col = "magenta")
axis(1,at=seq(1975,2015,5), labels=seq(1975,2015,5))
axis(2,at=seq(-40,80,20), labels=seq(-40,80,20))
box()
title(ylab = "US Dollars (constant) per capita\n2 year moving average",
      main = "Figure 6\nForeign capital flows to autocratic countries in SSA")
legend("bottomleft", inset = .01, c("Remitances","Aid","FDI"),lty=c(1,3,2), 
       col = c("red","black","blue"))


Year <- sort(unique(SSACONTROLLERS$year))
ProducerRemit <- NA
NonProducerRemit <- NA

for(i in 1:length(Year)) {
  Psubset <- SSACONTROLLERS[which(SSACONTROLLERS$year==Year[i]),c("ccode","remitraw","popraw")]
  Psubset$remit <- Psubset[,2]/Psubset[,3]
  ProducerRemit[i] <- weighted.mean(Psubset[,4],Psubset[,3], na.rm = TRUE)
  NPsubset <- SSANONCONTROLLERS[which(SSANONCONTROLLERS$year==Year[i]),c("ccode","remitraw","popraw")]
  NPsubset$remit <- NPsubset[,2]/NPsubset[,3]
  NonProducerRemit[i] <- weighted.mean(NPsubset[,4],NPsubset[,3], na.rm = TRUE)
}
ProducerRemitAVG <- SMA(ProducerRemit,5)
NonProducerRemitAVG <- SMA(NonProducerRemit,5)
#cbind(Year, ProducerRemit, NonProducerRemit)

plot.new()
plot.window(xlim=c(1975,2015),ylim=c(0,60))
lines(Year, ProducerRemitAVG, col = "red")
lines(Year, NonProducerRemitAVG, lty = 3)
axis(1,at=seq(1975,2015,5), labels=seq(1975,2015,5))
axis(2,at=seq(0,60,10), labels=seq(0,60,10))
box()
title(ylab = "US dollars per capita (5 year moving average)",
      main = "Figure 4.4\nRemittances in Sub-Saharan Afica")
legend("topleft", inset = .05, c("Currency Controllers","Non-controllers"),lty=c(1,3), col = c("red","black"))


Model <- plm(seigniorage ~ remit + gdppc + cbi,
             data = SSA, index=c("ccode","year"), model="random")
summary(Model)

stargazer(Model, style="apsr", no.space=TRUE, 
          single.row=FALSE,
          model.numbers = TRUE, title = "Table 4.4: Seigniorage Revenue", dep.var.labels.include = FALSE)

Failure <- glm(fail ~ cbi + remit + cbi*remit + party + gdppc + gdpgrowth + 
                   conflict + aidlog + oil + manufacturing + time*party +
                    time + time2 + time3 + trend + trend2 + 
                    factor(ccode), data = SSA, family = binomial(link = "cloglog"))

summary(Failure)
stargazer(Failure, style="apsr", 
          single.row=FALSE, model.numbers = FALSE, column.labels = "(7)", omit = c("time","trend"),
          notes = "cubic time polynomials and trend term not displayed",
          title = "Table 4.5: SSA Regime Failure", dep.var.labels.include = FALSE)

interplot(Failure, var1 = "remit", var2 = "cbi", xmin = 0, xmax = 1) +
  xlab("Central Bank Independence") +
  ylab("Estimated Coefficient on Remittances\nin survival analysis (logit regression)") +
  ggtitle("Figure 7\nConditional Effect of Remittances on Regime Failure") +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed")




Failure2 <- clogit(fail ~ remit + cbi + remit*cbi + party + gdppc + gdpgrowth + 
                     conflict + aidlog + oil + manufacturing + time*party +
                     time + time2 + time3 + trend + trend2 + 
                     strata(ccode), data = SSA)
summary(Failure2)

Failure3 <- glm(fail ~ cbi + remit + cbi*remit + party + gdppc + gdpgrowth + 
                 conflict + aidlog +  manufacturing + time*party +
                 time + time2 + time3 + trend + trend2 + 
                 factor(ccode), data = SSA, family = binomial(link = "cloglog"))

summary(Failure3)
interplot(Failure3, var1 = "remit", var2 = "cbi", xmin = 0, xmax = 1) +
  xlab("Central Bank Independence") +
  ylab("Estimated Coefficient on Remittances\nin survival analysis (logit regression)") +
  ggtitle("Figure 7\nConditional Effect of Remittances on Regime Failure") +
  theme(plot.title = element_text(hjust = 0.5)) +
  geom_hline(yintercept = 0, linetype = "dashed")


